#!/usr/bin/env python
from dolfin import *
from mshr import *
import sys


input_line = sys.stdin.read()
input_line.strip()
[par1, par2] = input_line.split()
par1 = float(par1)
par2 = float(par2)

x0=par1
y0=par2

resolution = 300

domain =   Rectangle(Point(0.0, 0.0), Point(1.0, 1.0)) \
         - Rectangle(Point(0.0, 0.0), Point(0.01, 0.3)) \
         - Rectangle(Point(0.99, 0.7), Point(1.0, 1.0)) \

domain.set_subdomain(1, Rectangle(Point(x0, y0), dolfin.Point(x0+0.1, y0+0.1)))

mesh = generate_mesh(domain, resolution)
domains = dolfin.MeshFunction("size_t", mesh, 2, mesh.domains())
#plot(domains, "Subdomains")
#interactive()


# Create classes for defining parts of the boundaries and the interior
# of the domain
class Left(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.01) and between(x[1], (0.0, 0.3)))

class Right(SubDomain):
    def inside(self, x, on_boundary):
        return (near(x[0], 0.99) and between(x[1], (0.7, 1.0)))

left = Left()
right = Right()

# Initialize mesh function for boundary domains
boundaries = FacetFunction("size_t", mesh)
boundaries.set_all(0)
left.mark(boundaries, 1)
right.mark(boundaries, 2)

# Define input data
a0 = Constant(1.0)
a1 = Constant(10000.0)
voltage = 1.0

# Define function space and basis functions
V = FunctionSpace(mesh, "CG", 1)
u = TrialFunction(V)
v = TestFunction(V)

# Define Dirichlet boundary conditions at top and bottom boundaries
bcs = [DirichletBC(V, voltage, boundaries, 2),
       DirichletBC(V, 0.0, boundaries, 1)]

# Define new measures associated with the interior domains and
# exterior boundaries
dx = Measure("dx")[domains]
ds = Measure("ds")[boundaries]
f = Constant(0.0)

# Define variational form
a = inner(a0*grad(u), grad(v))*dx(0) + inner(a1*grad(u), grad(v))*dx(1)
L = f*v*dx()

# Solve problem
u = Function(V)

# Define goal functional (quantity of interest)
#M = u*dx()
#tol = 1.e-4
#problem = LinearVariationalProblem(a, L, u, bcs)
#solver = AdaptiveLinearVariationalSolver(problem, M)
#solver.parameters["error_control"]["dual_variational_solver"]["linear_solver"] = "cg"
#solver.solve(tol)

solve(a == L, u, bcs)

# R = U/I
if(0):
	n = FacetNormal(mesh)
	m1 = a0*dot(grad(u), n)*ds(1)
	I1 = assemble(m1)
	print I1
	n = FacetNormal(mesh)
	m2 = a0*dot(grad(u), n)*ds(2)
	I2 = assemble(m2)
	print I2
	current = (I2+(-I1))/2
	print current
	R = voltage/current

# R = U*U/P
if(1):
	energy = a0*inner(grad(u), grad(u))*dx(0)
	E0 = assemble(energy)
	#print "Power 0 :", E0
	energy = a1*inner(grad(u), grad(u))*dx(1)
	E1 = assemble(energy)
	#print "Power 1:", E1
	power = E0+E1
	#print "Power total:", power
	R = voltage*voltage/power 

print R


# Plot solution and gradient
#plot(u, title="u")
#plot(grad(u), title="Projected grad(u)")
#interactive()
